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We present a new efficient method for computing the non-linearity parameters of the higher order 
correlation functions of local type curvature perturbations in inflation models having a A/"-componcnt 
scalar field, focusing on the non-Gaussianity generated during the evolution on super-horizon scales. 
In contrast to the naive expectation that the number of operations necessary to compute the n- 
point functions is proportional to N n , it grows only linearly in N in our formalism. Hence, our 
formalism is particularly powerful for the inflation models composed of a multi-component scalar 
field, including the models in which the slow-roll conditions are violated after the horizon crossing 
time. Explicit formulas obtained by applying our method are provided for n = 2, 3, 4 and 5, which 
correspond to power-, bi-, tri- and guarf-spectra, respectively. We also discuss how many parameters 
we need to parameterize the amplitude and the shape of the higher order correlation functions of 
local type. 

PACS numbers: 



I. INTRODUCTION 



Current observations of the cosmic microwave background (CMB) anisotropics indicate that primordial curvature 
perturbations are almost Gaussian [1J. In general, if the perturbations are purely Gaussian, the statistical properties 
of the perturbations can be completely described by the two-point correlation function (=power spectrum). On the 
other hand, if the perturbations deviate from the Gaussian distribution, the non-Gaussianity affects the higher order 
correlation functions, or higher order spectra. Currently, the non-Gaussianity is attracting attention as a powerful 
probe to discriminate various inflation models 0, I n particular, there are a large number of studies on the three- 
point correlation function (=bi-spectrum). However, the four-point correlation function (=tri-spectrum) can also 
be constrained by future accurate measurements [3, M, B Q- Using the analysis of both the bi-spectrum and the 
tri-spectrum in the future experiments, it is expected that we can extract more information about the mechanism of 
generating the primordial curvature perturbations. Hence, it is important to obtain useful formulas for the higher 
order correlation functions of primordial curvature perturbations. 

Roughly speaking, the leading order of the connected part of n-point function is 0(P n ^ 1 ), where P ~ 10~ 10 is the 
amplitude of the power spectrum. Hence, it is naively expected to be difficult to measure higher order correlation 
functions. However, when the non-Gaussianity is large, this estimate 0(P ra_1 ) will be replaced with 0(f^£ 2 P n ^ 1 ) 
or even larger. Here, /nl is a non-linearity parameter given in Refs. @, E|; 

C = Cg + |/nlCg, (I-l) 

where £ is the curvature perturbation on uniform energy density hypersurface and (g is the linear Gaussian part. No- 
tice that observationally /nl can be as large as 0(100). This possible enhancement slightly improves the detectability 
of the higher order correlation functions. Furthermore, the number of argument wavenumbers of the n-point function 
is n — 1. When the CMB temperature anisotropies, Ci, are measurable up to I — £ max , the number of independent 
wavenumbers which we can measure will be roughly estimated as ^, 2 nax . Hence, the number of different combinations of 
argument wavenumbers increases as imSx ■ This large number enhances the effective amplitude of n-point function 
to 0(f^(fj<sh^maxP) n ~ 1 ), while the amplitude of Gaussian noise is 0{P n / 2 ). The detectability of the n-point function 



* Electronic address: 
i Electronic address: 

* Electronic address: 



shu@a. phys.nagoya-u.ac.jp 
suyamaOicrr. u-tokyo .ac.jp 
tama@scphys.kyoto-u.ac.jp 



2 



is basically determined by the ratio of these two numbers, O (f N ^P (f'NL^maxVP) n 1 j- Hence, if f^L^maxVP 
exceeds unity, all the higher order correlation functions are in principle measurable. For the Planck satellite Q, it 
is expected that £ m ax ~ 0(2000). Hence, naively, if /nl would be as large as O(50), f^imax^/P can exceed unity. 
This fact strongly motivates a systematic derivation of the formulas for higher order correlation functions. 

In this paper, we present a new method to calculate general rt-point functions of local type primordial curvature 
perturbations. This new method is much more efficient than the straightforward calculations, especially when applied 
to the models with many components of inflaton field, including the models in which the slow-roll conditions are 
violated after the horizon crossing time. This method is based on the diagrammatic approach given in Ref. Q 
as well as on our previous work [Tol . [Tl| . in which the formulation for the bi-spectrum was developed. As for the 
parameterization of the higher order spectra, it is well known that the bi-spectrum can be parameterized by a single 
parameter, so-called non-linearity parameter, /nl, while the tri-spectrum is parameterized by two parameters tnl 
and gML |12l ] due to the existence of two distinct terms that exhibit a different wavenumbcr dependence. That is, 
the number of parameters necessary to describe the higher order correlation functions is equal to the number of 
independent terms which have a different wavenumber dependence. Based on the diagrammatic method, we also 
show that one can easily count how many parameters we need to parameterize the amplitude and the shape of higher 
order spectra of local type. 

This paper is organized as follows. In section [II] we briefly review the SN formalism [l3|, [lj, [ID, [HI, [TtJ , which is the 
foundation of our present analysis. We also discuss how many parameters we need to parameterize the higher order 
correlation functions. In section HTT1 we present our diagrammatic method for the computation of n-point correlation 
functions of primordial curvature perturbations. As an application of our method, in the succeeding section IIVI we 
give concise formulas for the power-, bi-, tri- and quad-spectra, of the primordial curvature perturbations generated in 
multi-component inflation models. Section [V] is devoted to discussion and conclusion. 



II. LOCAL TYPE PRIMORDIAL CURVATURE PERTURBATIONS AND THEIR 

PARAMETERIZATION 

We focus on the non-Gaussianity generated during the evolution on super-horizon scales in multi-scalar inflation. 
We start with a brief review of the SN formalism. Using the SN formalism, we present a diagrammatic representation 
for general n-point functions of local type primordial curvature perturbations, and show how they are parameterized. 

A. Background equations 

We consider a A/"-component scalar field whose action is given by 



S = - I d 4 x^ 



hiug^d^d^ + V{ct>) , (11.1) 
(J,J=1,2,--- ,A0 , 

where g^ is the spacetime metric and hu is the metric on the scalar field space. In this paper we restrict our 
discussion to the flat field space metric hu = 5j,j to avoid inessential complexities due to non-flat field space metric, 
though the generalization is straightforward 
We define ip\ (i — 1, 2) as 1 

<p[ = 4, 1 , V \ = 4>\ (II.2) 

where a dot " ' " represents differentiation with respect to the cosmological time. 

For brevity, hereinafter, we use Latin indices at the beginning of Latin alphabet, a, b or c, instead of the double 
indices, i.e., X a = X\. Then, the background equation of motion for ip a is 

^ a = F a { V ) , (II.3) 



Here, we take different definition for i/^ from that introduced in our previous papers [lot ITU , which was defined as </?2 = defr 1 /dN. 
Based on previous definition, specific expressions for P a b or Qf^^ ...{, 1 defined as Eq. JIII.3I in the later Sec. IHI Xl include the terms 

which diverge when V = 0, which is not a suitable formulation for the numerical calculations. If we define y>\ as in Eq. QL2) , there are 
no divergences of P a . or Q?„,, , . at the time when V = 0. 
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where TV is the e-folding number and F a (= F[) is given by 



with = 6 IJ (dV/d(f> J ). The homogeneous background Friedmann equation is given by 

H 2 = \ (\<pi<P2i + V 



(II.4) 



(II.5) 



with ip 2 i = Sij<P2- 

In the SN formalism [I3, LLJ, [la, Us j the difference in e-folding number between two adjacent background 
solutions describes the evolution of £, curvature perturbations, on super-horizon scales. The solution of the background 
inflationary dynamics dominated by a TV-component scalar field is labelled by 27V— 1 integration constants A a , besides 
the trivial time translation 5N. Let us define 6ip a as the perturbation, 



<fy>°(A; TV) = V a (\ + SX; N) - V a (\; TV) , 



(II.6) 



where A is abbreviation of A a and <5A a is a small quantity of 0(5). 

27V parameters {TV, A a } parameterize the initial values of fields. There is an arbitrariness in choosing the integral 
constants, i.e. a different choice of integration constants A a = / a (A) is equally good. Here we leave the choice of A a 
unspecified since all the discussion in the paper is not affected by the choice. 

5ip a (TV) defined by Eq. (|II.6[) represent perturbations of the scalar field on the TV = constant gauge . In this 
case, C at each point in space depends on the fluctuations of the scalar field at the same spatial point, and is given by 



C(iV F ,5?) = ^iA, 



(II.7) 



TV* 

ai_a-2... a 7 



d n N(N F ,{<p a }) 



dtp ai dip a2 ■ ■ ■ dtp a ™ 



(0) 



Here, the values of scalar fields on the initial flat hypersurface, {<p%}, differ from place to place and characterize the 
initial perturbation. Since the e-folding number between the initial flat hypersurface and the final uniform energy 
density hypersurface depends on {</?"}, as its argument we have used {ip a } instead of the initial time TV = TV*. We 

(0) 

have decomposed the scalar field as ip a = f a + 5ip a and Taylor expanded C in terms of 5ip a = 0(6). The suffix 
* represents the value evaluated at a certain time TV* which is shortly after the horizon crossing time. The final 
hypersurface at TV = TVf is chosen to be an uniform energy density surface. As is well known, C(7Vf) is independent 
of the choice of TVf as long as TVf > TV C , where TV C is a certain time after the background trajectories have completely 
converged. According to the 5N formalism, the expansion coefficients TV* ia2 a ^ are simply given by the derivatives of 
N(Np, {ip a }), where N(Np, {<p a }) is the e-folding number spent during the evolution of the homogeneous universe, 
in phase space, from the initial point {ip a } to the final uniform energy density surface. 



B. Parameterization of the n-point functions 

Let us begin with the two-point function. At the leading order in S, 2 
the two-point function of £ can be written as 

(Cg 1 Cfe>o = Pdki)S (3) (ki + k2) 

= NaNbA^iki^Hki + ha) , (H.8) 

where (• • • ) c means the expectation value of the connected part of "• • • ", and we have abbreviated the suffix *. Here 
we have introduced the covariance matrix A ab (k) defined by A ab (k)6(k + k') = (5ip a j.Sip b -,) c . 



2 In Rcfs. 9, 18], the authors have also considered the one-loop corrections which are the higher order in 5. In this paper we consider 
the only tree-level spectrum and neglect the one-loop corrections. 
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We assume that all the relevant components of the scalar field satisfy the slow- roll conditions at least until N = , 
in which our formalism works quite efficiently. Otherwise, correlation functions can not be parameterized by a small 
number of parameters. In this case, {Sip d } is approximated by a set of Gaussian random variables with the scale 
invariant spectrum 3 , and A ab is given by A ab — A ab P(k) with 

A[( = 5" , A{i = AH = A" = , P(k) = ^ 2 . (II.9) 

Strictly speaking, even in the slow-roll inflation, 5(p% deviates from pure Gaussian perturbation due to the effect of 
interaction. However, the non-Gaussianity of C caused by this deviation is suppressed by the slow-roll parameters, 
which is an undetectable level in the future experiments pjj H3, HH, [22J. Hence, we neglect the non-Gaussianity of 
S(f>l here. 4 In a similar fashion, the three-point correlation function at the leading order is given as [13] 

(C&CfeCg^c = 777— # (Pdh)Pc(k 2 ) + 2 perms.) S^(h + k 2 + fc 3 ) , (11.10) 

where A a = A ab Ab. In deriving this equation, we have used {5ip a t 5tp b - Sep * ) c = 0. From this equation, we find 

that (C^C^Cfcglc is 0(5 ). Since the wavenumber dependence of the bi-spectrum is completely given by the products 

of the power spectrum, the bi-spectrum is characterized by a single parameter N a N b N ab /(N c N c ) 2 , which controls 
the overall amplitude. Following literatures, we redefine the non-linearity parameter /nl given in the introduction 
as [13 

S N a N b N y. 

/nl = ^ V JV f . (11.11) 

If only one field contributes to the curvature perturbation, /nl defined by Eq. (jll.lip is equivalent to Eq. 1|) . While 
Eq. (II. lj) is valid only when the single field dominates the curvature perturbation, Eq. (jll.lip can be applied to larger 
classes of inflation models where the curvature perturbations are sourced by multiple fields. 

We can also write down the leading order four-point correlation function (the tri-spectrum) as fl2l . Il7j 



N a N ab N bc N c ( p c(fci) p c(fci2) p c(fc4) + n perms } 
(N d N d ) 



N a N b N c Nabc 

(N d N d f 



(P c (fc 1 )F c (fc 2 )P c (fc 3 ) + 3 perms.) 



^ 3 )(fci + fc 2 + fc3 + fc 4 ) , (H.12) 



where kij = \ki + kj\ and N ab = A ac A bd N c d- We see that the four-point function is 0(8 e ). Unlike the bi-spectrum, 
the tri-spectrum has two distinct terms that exhibit different wavenumber dependence. As a consequence, we need 
two parameters to specify the tri-spectrum. Following Ref. [T^ |. we use the non-linearity parameters tnl and (?nl 
defined by 

N a N ab N bc N c 25 N a N b N c N abc 

(N d N d f 54 (N d N d ) 3 

We can further proceed to higher order correlation functions. The main issue that we address in the rest of this 
section is how many parameters are necessary to parameterize the local type n-point function. Of course, we can 
count the number of such parameters by directly calculating the n-point function from Eq. (III.7j) as in the case of the 
bi-spectrum or the tri-spectrum. Although in principle there are no difficulties in such a direct counting, the actual 
computation becomes exponentially more cumbersome as we proceed to higher order. Here, instead of resorting to 
the direct computation, we use a diagrammatic method 



3 As is well known, the deviation from the scale invariant spectrum can be given by the slow-roll parameters at the horizon crossing 
time of the corresponding scale, and , A^f and A22 are also suppressed by the same slow-roll parameters. Since we know that 
the deviation from the scale invariance is observationally small, it is natural to assume that the slow-roll conditions are well satisfied 
at around the horizon crossing time. Therefore, for simplicity, we evaluated A ab at the horizon crossing time in the limit of vanishing 
slow-roll parameters. 

4 Note that in the slow-roll inflation the non-Gaussianity of <5</?J dominates the trispectrum of f . But it is too small to be detectable in 
the future experiments 0, HO, H3, HI ■ 
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The leading order of the n-point function consists of terms of 0(5 2n ~ 2 ), which are given by products of (n — 1) 
power spectra. According to the diagrammatic method, each of these leading terms has a corresponding connected 
diagram that consists of n vertices and (n — 1) lines connecting two vertices. Such a connected diagram should have 
a tree structure. Namely, there is always a unique path that connects any pair of vertices in the diagram. We refer 
to such diagrams as reduced tree diagrams, to distinguish them from the (full) tree diagrams that will be introduced 
later. 

The rules of the reconstruction of the leading term that constitutes the n-point function from a given reduced tree 
diagram are as follows [9j. First, we assign a different wavenumber ki(l < i < n) to each vertex • of the diagram, 
where ki, k 2l ■ ■ ■ ,k n are the arguments of the n-point function with the constraint ki + k 2 + ■ ■ ■ + k n = 0. Next 
we assign a wavenumber to each line in the diagram, too. In general, removing a line from the diagram yields two 
respectively connected sub-diagrams. Then, one assign to the removed line the sum of the vectors associated with all 
vertices in one of the two sub-diagrams. We do not care which of two sub-diagrams we choose since only the length 
of the assigned waved number is used in the following discussion. An example of the assignment of the wavenumber 
is given in Fig. [TJ 




FIG. 1: This figure shows the reduced tree diagram corresponding to one of the leading terms which constitute the n-point 
function whose arguments are fci, fez, • • • , k n - The assignment of wavenumbers to the vertices and lines is illustrated. 

After associating the wavenumbers with all lines, now we can assign the corresponding factors to the vertices and 
the lines. As for the vertex with p lines attached, assign the factor N aia2 ... a to it. As for the lines, assign A ab P, 
where the argument of the power spectrum P is set to the length of the wavenumber associated with each line. By 
multiplying all these factors assigned to vertices and lines, and summing up all independent diagrams which are not 
mutually isomorphic, wc obtain a function of n wavenumbers, which constitutes the n-point function. The indices 
in N aia2 ... ap assigned to each vertex are contracted with the indices of p neighboring lines. Contraction is performed 
between lower and upper indices as usual. 

Here, we did not associate a factor l/p\ with the vertex N aia2 ... ap from the beginning for the following reason. A 
vertex with p lines attached has p lower indices to be contracted with the upper indices in A ab P associated with the 
p attached lines. These p lines are all to be distinguished because they are all labelled with different wavenumbers. 
Therefore there are p\ ways of contraction between two sets of p indices. If we do not distinguish which indices are 
contracted, the factor 1/pl associated with the vertex is canceled. 

Finally, by taking the sum over all the possible reduced tree diagrams, we obtain the n-point function. As an 
illustration, we show the diagrams for n = 3 and 4 in Figs. [2] and [3l respectively. 

It is not a trivial matter whether the functions constructed from two reduced tree diagrams that are not isomorphic 
to each other always yield a different functional dependence on the wavenumbers. As we explained in the appendix [Al 
if the two reduced tree diagrams with n vertices are not mutually isomorphic, the corresponding functions of n 
wavenumbers are always different [23[. Therefore, the number of parameters necessary to determine the n-point 
function of £ is equal to the number of independent reduced tree diagrams with n vertices. As an illustration, we show 
in Table [J the number of free parameters and the corresponding diagrams for n = 3, 4, 5 and 6. 

A similar diagrammatic approach for the higher order correlation functions has been developed in the context 
of galaxy correlation or large scale structure. In Ref. [24[, the author has given a numeration of the number of 
independent tree diagrams for general n, using the generating functions based on the combinatorial analysis [25|. 
Applying this method to our discussion about the higher order correlation functions of the primordial curvature 
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A ab P{k 1 ) 



N d 



A cd p(k 3 ) h 



FIG. 2: This diagram represents the leading term of the bi-spectrum of primordial curvature perturbations. 



n = 4 



k 



I A be P{k 2 ) 



N a 

\ 

A od P(fei) 
ku 



k. 



N. 



def 



A cf P{h) at 



A ab P(k 1 

K 

N bc 



N f 

k 4 e> 

A ef P(k 4 ) 

k 3 

A'' d P(k l2 ) % de 



FIG. 3: These diagrams represent the leading order terms of the tri-spectrum. These two distinct diagrams show different 
wavenumber dependence. Hence, two parameters are needed to describe the tri-spectrum. 



perturbations, we can find the number of independent reduced tree diagrams for general n, which corresponds to the 
number of free parameters for general n-point functions. 



III. A NEW METHOD TO COMPUTE n-POINT FUNCTIONS 



In this section, we will provide an efficient method to compute the non-linearity parameters to characterize the 
n-point correlation functions. 



A. tree-shaped diagrams 

We start with the fact that C(-^f) is independent of the choice of the time of the initial flat hypersurface iV*. By 
choosing N* to be identical to Np in (|II.7|) . we obtain [l6[ 

C(^) = E^<a,..a„^^ 2 ---^ , (HIT) 

where Sip F — 5ip a (Np) are field perturbations evaluated on the flat slice at N = Np. As we have mentioned in the 
previous section, C(-^Vf) represents the curvature perturbation on the uniform energy density slice at N — Nf. Hence 
the above equation (jlll.ip means that the curvature perturbation £(Np) is simply caused by a time shift between the 
flat slice and the uniform energy density slice at the final time TV = Np. Hence, as shown in the appendix[Bl N^ iCl2 a 
can be written only by local quantities at N = Np. N^ a2 a are therefore obtained immediately, once we specify 
■■P'l- 
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number of 
free parameters 


corresponding reduced \xee diagrams 
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1 (fm) 
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4 
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I . . I. I I . v-l 



TABLE I: This table shows the number of free parameters and the corresponding reduced tree diagrams for n-point correlation 
functions with n = 3, 4, 5, 6. 



What we need to evaluate is {8p F } as functions of {Sip®}. The evolution equations for Sp a , which can be obtained 
by perturbing the background equation (|II.3[) . are given by 



^S<p a (N) = P a b Sp b (N) + ±Q? 3)bc (N)Sp b (N)Sp c (N) + 



1 



where P a b and Q?. 



are, respectively, defined by 



dF a 

pa _ " x 



(0) 

<p=<p(N) 



gi-1 pa 



dip bl dip b2 ■ ■ ■ dip be - 



(0) 

<p=tp(N) 



(111.2) 



(III.3) 



For the purpose of the evaluation of the n-point function, it is enough to truncate the expansion on the right hand 
side in Eq. (|III.2[) at (n — l)-th order. By solving the above equations from N = N* to Np with initial conditions 
Sip a (N*) = Sipl, we obtain Sip F expressed in terms of {Sip%}. Due to the non-linear evolution after the horizon crossing, 
the distribution of {Sp F } is in general non-Gaussian even if that of {Sip*} is Gaussian. 

If we solve the equations (|III.2I) iteratively, we can express formally Sp F as a Taylor expansion in terms of {5<p%}. Let 

us denote the m-th order terms in the iterative expansion by 5 V p. Then Sp a F can be written as 6ip F — X)m=\ 
where we truncate the expansion at the (n — l)-th order because higher order terms are irrelevant to the 7i-point 

function. By definition, S^P a F contains m Gaussian random variables, {Sip"}. Namely, there is a factor Sip" 1 ■ ■ -Sp" m 



in Sfp. The indices in this factor are to be contracted with the interaction vertices Q?t)b 1 b x ...bi 
function A a b , which obeys 



or the Green 



— A%(N,N') = P a c (N)A c b (N,N>) , 
^La%(N,N>) = -A" c (iV,iV')P c h (iV') 



(IIL4) 
(III.5) 



Here again an upper index is contracted with a lower index, as usual. All the possible ways of contraction contribute 

S'a. 



to S%\. 



We can associate a diagram as presented in Fig. @] with each way of contraction. Hereinafter, we refer to such a 
diagram as a tree-shaped diagram, to distinguish it from the reduced tree diagram introduced earlier and from what 
is simply called a tree diagram which will be introduced later. A tree-shaped diagram is drawn obeying the following 
simple rules. We start with a solid circle • and attach a line downward to it. We attach an interaction vertex ® to 
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V ~ 
The total number of Sip" is m 

(m < n — 1) 



FIG. 4: One example of tree-shaped diagram, which contributes to 8 V 5 The diagram branches off from the top to the bottom. 



the other end of this line. From the vertex, several lines extend downward and they end with another interaction 
vertex ® or a half open circle. This process is repeated until all the end points are terminated by a half open circle. 
The total number of half open circles should be m. 

The solid circle • corresponds to , and hence we assign the time N = Np to it. A half open circle corresponds to 
the initial Gaussian variables Sip®, and hence the time N = N* is assigned to it. An interaction vertex ® corresponds 
to the factor Q/£\ bl ... be where N (iV* < N < Np) is the time assigned to this vertex and £ (3 < I < m) is the 

number of attached lines. Here, we did not associate a factor l/(£ ~ 1)1 with Q b ^ bl ... b( 1 (N). The reason is the same 
as before in the case of N aia2 ... ap discussed in Sec. Ill Bl Here, the half open circles are all supposed to be labelled, i.e. 
distinguishable. In this case dropping the factor 1)! is exactly compensated by not distinguishing the order of 

lines attached to the same interaction vertex. 

In this diagram the time flows from the bottom to the top as is indicated in Fig. [4] Finally, to each line segment, we 
assign the Green function (propagator) A a b (N, N'), where N and N' are the time coordinates assigned to the upper 
and lower ends of the line, respectively. Contracting upper and lower indices between the adjacent objects (the solid 
circle, interaction vertices, half open circles and lines) and integrating over all the time coordinates assigned to the 

interaction vertices in the whole range of their possible variation, we obtain a quantity which constitutes N„ Sf F . 

Collecting all terms corresponding to different diagrams yields total N^S <P F . 



B. n-point correlation functions 

Instead of ((Np), we first compute n-point functions of Cp defined by 

Cf^A^, (III.6) 

which is the linear truncation of the Taylor expansion of £(Np) in terms of 5(p F . The n-point function of Cf'"^ is 
given by the sum of all the possible connected tree diagrams obtained by contracting all the half open circles {S<p®} 
in pair from a product of n tree-shaped diagrams. (See Fig. [5]) Contraction between the half open circles within the 
same tree-shaped diagram can be neglected because it produces a loop. For the same reason, there is not more than 
one contraction between any pair of tree-shaped diagrams. Let us represent this contraction between a pair of half 
open circles, by a full open circle o, to which A ab defined in Eq. (|II.9[) is assigned. We refer to the diagram obtained 

by this contraction simply as a tree diagram. The leading terms of the n-point function of £p are 0(S 2 ^ n ~ 1 " 1 ), and 
hence the tree diagram should have (n — 1) open circles o. As any pair of tree-shaped diagrams does not have more 
than one contraction between them, all the half open circles belonging to a single tree-shaped diagram are contracted 
with different tree-shaped diagrams. Since all the tree-shaped diagrams are labelled with a different wavenumber fcj, 
the assumption that all the half open circles are distinguishable holds. This contraction process does not produce any 
further statistical weight, and hence all the tree diagrams have the same weight of unity. 
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power spectrum bi-spectrum 

FIG. 5: Tree diagrams corresponding to the power- and bi-spectra. 

In addition to the linear term qJ , C(Nf) contains terms non-linear in 8ip F , which also contribute to the n-point 
functions of ((Np). There is one-to-one correspondence between the non-linear terms in ((Np) and the interaction 
vertex <g> that is directly connected to • by a line without any intervening vertices o or ®. Hence, we can take into 
account this non-linear contribution simply by replacing the interaction vertices directly connected to • as 

N[A a c (N F ,N)Q c ie)bi ... be i (N) N^A a c (N F ,N)Q c Wbi ... be _ i (N) 

= N^A\(N F ,N)Q c {i)bi ... b( _ i (N)+N b F i ... be _ i 5 (N — (Np — e)), (III.7) 

where e is an infinitesimally small number. By this prescription, we can obtain the non-linearity parameters only 
from the tree diagrams. 

Now we are ready to show that our method to evaluate the n-point functions can be reduced to the problem of 
solving the ordinary differential equations for vector variables that have only a single index a (1 < a < 27V), which 
is the main result of this paper. Let us consider one tree diagram which constitutes the n-point function. We focus 
on one of sub-diagrams obtained by removing one vertex ® or o from a tree diagram, which we denote by T a or 
r a . If the line which was attached to the removed object is pointing downward (upward), the vector has a lower (an 
upper) index. Suppose that the object attached to the other end of this line is an interaction vertex (g) with which 
Q(e)bi---b e iC^O associates. Let us consider the case of a vector T a with an upper index. Notice that the other lines 
connected to this vertex are also similar sub-diagrams which consist of smaller number of vertices than that we are 
focusing on. We denote the product of the vectors associated with all these sub-diagrams by a tensor M Cl ' " c< - 1 . Those 
vectors appearing in M 01 '" '- 1 are already known by the induction assumption. Then, T a can be defined recursively 
as 



IN, 

From this equation, we find that r a (7V) satisfies 



a (N) = f N dN> A a b (N, N')Q\ t)ci ... Ct i {N')M*-°'-i (N 1 ). (III.8) 



A-f a {N) = P%(N)f b (N)+Q a Wci ... Cei (N)M c --^(N). (III.9) 

The boundary conditions for r a (iV) are set by F a (A?*) = at N — N*, and hence the above equation is to be solved 
in the forward direction in time. 

In the case of the vector with an upper index r a (7V), the neighboring object can be o instead of ®. In this case 
the initial conditions are given by ^(TV*) = A ab Tb(N f ), where Tb(N r ) is the vector corresponding to the sub-diagram 
with the neighboring vertex o being removed. The equation to solve is simply the homogeneous one given by 

Ar a( iv) = -r b (N)P b a (N). (iii.io) 

Similarly, for a vector T a with the neighboring object being ®, we have 

fN F 

T a (N) = / dN> A-(/v',7V)g^ )ci ... c ^ i (iV')M 6 C2 - c '- 1 (iV') , (111.11) 

J N 



and this vector obeys 



d T a (N) = -T b (N)P b a (N) Q\ k)aC2 ... Ct _SN)M b c ^ c ^{N) . (111.12) 



dN 
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The boundary conditions for T a (N) are set by r a (7Vp) = at N = Np, or equivalently T a (NF — e) = N aC2 ,.. Ce _ 1 
M C2 "' Cl - 1 (Np) taking into account the <5-function term in the definition of Q%\ bl ... be i in Eq. (|III.7[) as boundary 
conditions. In this manner, the effect of the non- linear terms in C(-^f) in 1|) can be absorbed by the boundary 
conditions in general. The equation is solved backward in time. There is another case in which the neighboring object 
is •. This simplest case can be also handled in a similar manner. We defer its explanation to the succeeding section, 
where we exhibit some more explicit formulas. In Tabic HH we summarize the notation of the vector quantities which 
will be used below, showing the correspondence to the tree-shaped diagram. 



vector quantities 

r„(iV;diag), f°(W;diag) 


diag 


trees-shaped/ 
diagrams 


T„(N;A) = N a (N) 


A 


/ 


J Aa ( AT. ~A \ AT a ( AT'\ 
1 \pi \ J%) = iV [1\ ) 






T a (N;B) = n a (N) 


B 




r a (N;B) = Q"(N) 


B 




T„(N;C) = $„(JV) 


C 




T a (N;V) = * (AT) 


V 




r"{N;£) = fT'(iV) 


£ 





TABLE II: This table shows the summary of the correspondence between the vector quantities used in this paper. 



To obtain an expression for the n-point function written in terms of such vectors, we arbitrarily choose one vertex 
® or o from a tree diagram at the beginning. Suppose that the chosen vertex is an interaction vertex ® with which 
Qft)bi— b«_i(-^0 associates. After preparing all the necessary vectors, T a and T/L (1 < i < I — 1), which correspond 
to the sub-diagrams obtained when this vertex is removed, we can immediately write down the contribution to the 
n-point function from this tree diagram as 



N F 



dN r a (N)Q a (k)bi ... be _ i (N)]Jf b ( ; ) (N) . 



(III. 13) 



If the vertex which we initially focused on is o, with which A ab associates, we do not need the final integration over 
N. We denote the vectors that correspond to the sub-diagrams obtained by removing this o by ri 1} and r[ 2) . Then, 
we compute 

A ab r^(N,)T^{N,), (111.14) 



instead of the expression (|III.13[) . In this diagrammatic method, the final expression for the spectrum in appearance 
depends on which vertex we chose at the beginning, but, of course, all different looking expressions are equivalent. 
Practically, it is more efficient to choose a vertex near the center so as to reduce the number of necessary vectors, 
although the definition of the center of a diagram is not so clear in many cases. 



C. Relation to the reduced tree diagrams and statistical weight 

As we mentioned in Sec. Ill Bt the reduced tree diagram is useful to classify the wavenumber dependence of the 
n-point functions, while the (full) tree diagram is a powerful tool for explicit computation of the n-point functions. 
We show that the reduced tree diagram introduced in Sec. Ill Bl is actually a simplified version of the tree diagram. 
It will be manifest that the solid circles attached to the ends of diagrams have the same meaning in both diagrams. 
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In the reduced tree diagram the internal lines represent power spectrum of the initial Gaussian random field {Sip"}, 
which is expressed by an open circle o in the tree diagram. Hence, each line in the reduced tree diagram corresponds 
to a line with an open circle o in the tree diagram. The sub-structure described by the interaction vertices (£> in the 
tree diagram is completely abbreviated in the reduced tree diagram. Hence, there is a degeneracy such that different 
tree diagrams contribute to the same reduced tree diagram. As explained in Sec. Ill Bl and proven in appendix [XJ the 
wavenumber dependence of the n-point functions is classified by the topology of the reduced tree diagram. This means 
that plural tree diagrams can give the contribution to n-point function with the same wavenumber dependence. 

In Fig. [51 as an example, we show the diagrams corresponding to the tri-spectrum coefficient gNL- We can decompose 
the top-left reduced tree diagram into 4 sub-diagrams by cutting all lines off. These sub-diagrams are counter parts 
of the tree- shaped diagrams. The lower part of this figure explains correspondence between these sub-diagrams and 
the tree-shaped diagrams. There are two tree-shaped diagrams with four half open circles as is explicitly shown in this 
figure. Hence, we find that the formula for g^L is composed of two different terms. 

When we consider the statistical weight of the diagram, this correspondence between the reduced and full tree 
diagrams is important. The starting point is the fact that the statistical weight of each tree diagram is unity when 
each end point • is labelled by the assigned momentum. Therefore counting the statistical weight by writing down all 
different tree diagrams is straightforward. However, the non-linearity parameters are defined based on the reduced tree 
diagram. The most of the patterns which occur as a result of permutation of the momenta assigned to the end points 
• is taken care already in the definition of the non-linearity parameters. (See Eqs. fill. 10[) and (|II.12[> .) However, 
here we should notice that some of the half open circles in the tree-shaped diagrams can be distinguishable, while 
the sub-diagrams obtained from the reduced tree diagram as mentioned above do not distinguish their legs at all. 
Therefore when there are several distinguishable patterns to assign the labels to the half open circles in a tree-shaped 
diagram, the term containing such a tree-shaped diagram has a factor corresponding to the number of patterns. 

As an example, we again consider the case of gNL- The tree-shaped diagram that has two three-point interaction 
vertices shown in Fig. [6] has three distinguishable patterns in assigning the labels {1,2,3} to the three half open 
circles. This means that we need to add the corresponding factor 3 to the contribution containing this tree-shaped 
diagram. (See the expression in Eq. pV.6[) below.) 



9NL decomposing X 3 



reduced- diagram 



- 1 



correspondence 

< n — 




"decomposed" 
diagram 



and 



tree-shaped/ diagram 




FIG. 6: This figure shows the relation between the tree-shaped diagram and the reduced diagram for the non-linearity parameter 
gNL- The upper one shows how to decompose the reduced diagram, and the lower one shows the correspondence between the 
"decomposed" diagrams and the tree-shaped diagrams. 

We want to emphasize that the above formulation simplifies the computation of higher order correlation functions a 
lot. In this formulation, we have only to solve vector quantities with only one index. Therefore our computation scheme 
requires the number of operations proportional to Af in computing the non-linearity parameters in n-point function. If 
we performed a naive straightforward calculation, in which the derivatives of the e-folding number N are computed by 
using the finite difference method numerically, the required number of operations is proportional to Af n . If we naively 
performed perturbative expansion, in which we connect the interaction vertices by propagators A a b (7V, N') and perform 
integration over the time coordinates of the interaction vertices, the necessary number of operations would be even 
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7V„(7V) A 

FIG. 7: This figure shows how we decompose the tree diagram for the power spectrum, when we focus on the vertex o. 

larger. When the number of inserted interaction vertices is m, naively we need (3(several x (number of time steps)™) 
operations to compute the contribution of the single diagram. On the other hand, in our formulation we need only 
0(several x (number of time steps) x (to + a few)) operations. Therefore our scheme is particularly useful for the 
computation of higher order correlation functions in the inflation models with a large number of field components. 

IV. EXAMPLES 

In this section we apply our formalism to the computation of the power-, bi-, tri- and gwad-spectrum to demonstrate 
the efficiency of our method. 



A. Power spectrum 

Let us first consider the power spectrum. There is only one tree diagram that contributes to the power spectrum, 
which is shown on the left hand side in Fig. [5] Following the prescription given in the previous section, we focus on 
a unique vertex o. 

Then, we can decompose this tree diagram into this open circle, with which A ab associates, and two identical sub- 
diagrams shown on the right hand side in Fig. [7] Corresponding to this simplest sub-diagram, we introduce a vector 
N a (N), whose explanation was deferred in the preceding section. This vector is defined by the equation 

-^N a (N) = -P b a (N)N b (N), (IV.l) 
with the boundary conditions N a (Np) = . The vector N a (N) represents the derivatives of the e- folding number 

(0) I 1 

with respect to <p a evaluated at ip a = ( P a (N). Using this vector, the power spectrum of ((Nf) is expressed as [llj 

% = A ab N a {N*)N b (N*) = W m , (IV.2) 
where P( and P are those which have already appeared in Eq. III. 8|) . 



B. Bi-spectrum 



Just like the power spectrum, there is only one tree diagram that contributes to the bi-spectrum, which is presented 
on the right hand side in Fig. El Let us focus on the interaction vertex ® to which Q? 3 )t, c (N) is assigned, and 
decompose the diagram into the chosen vertex ®, a sub-diagram denoted by N a (N) and the two same sub-diagrams 
denoted by N a (N) as illustrated in Fig. M 

The sub-diagram denoted by N a (N) is reduced to that denoted by N a (N) if we remove one open circle o. Hence, 
following the general rule explained in the preceding section, the new vector N a (N) is obtained by integrating 



^N a (N) = P%(N)N b (N) , 



(IV.3) 



from N = N* with the initial conditions N a (N*) = A ab Nb(N*). Applying the general formula (|III.13j) supplemented 
by (|III.7[) . we recover the result previously obtained in Ref. [llj . 



(i 



NL 



w: 



N[ b N a (N F )N b (N F ) 



dNN a (N)Q a {3)bc (N)N b (N)N c (N) 



(IV.4) 
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T 



N a (N) 




x2 

N"(N) 



FIG. 8: This figure shows how we decompose the tree diagram for the bi-spectrum, when we focus on the interaction vertex i 



C. Tri-spectrum 



As we mentioned in the previous subsection IIIB| we need two parameters, tnl and gNL, for the tri-spectrum. 
The tree diagrams for the tri-spectrum were shown in Fig.[5J From this figure, we find that gNL consists of two tree 
diagrams. 



four-point function (tri-spectrum) 




reduced^ 
diagrams 



FIG. 9: The diagrams for the tri-spectrum. As we have mentioned in the previous section III Bl for the tri-spectrum due to 
the difference of the scale dependence we need the two parameters tnl and gNL, which can be also distinguished by using the 
reduced diagram. Using the tree diagram gNL is decomposed to two diagrams as shown in this figure. The arrows indicate the 
focused vertex in each diagram. 

We choose a focused vertex in each diagram as indicated by arrows in Fig. [5] Following the prescription explained 
in the preceding section, we have 



tnl = wc 3 [A ab n a (N*)n b (N*)] 
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g NL = ^W- 6 \ N* bc N a (N F )N b (N F )N c (N F ) + / dNN a (N)Q a {i)bcd (N)N b (N)N c (N)N a (N) 



N F 



3 j dNn a (N)Q a (3)bc (N)N b (N)N c (N) 



where Q a (N) is a new vector obtained by solving 

d 



dN 



n a (N) = -n b (N)P b a (N) - N b (N)Q b 3)ac (N)N c (N) , 



(IV.5) 



(IV.6) 



(IV.7) 



backward in time from N = N F with the boundary conditions Q a (N F ) = N^ b N b (N F ). The first (second) line on the 
right hand side of <?nl represents the contribution of the left (right) tree diagram corresponding to gNL in Fig- El In 
our formulation, it is enough to solve differential equations for three vectors N a , N a , fl a to compute tnl and <?nl- 
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D. guaci-spectrum 

In order to demonstrate the efficiency of our formulation, we show the explicit formula for the gtiad-spectrum. 
We also show the correspondence between the reduced diagram and the tree diagram for the fifth-order spectrum 
(five-point correlation function) in Fig. 1101 Using the formula, we obtain an expression for the gwarf-spectrum as 

five-point function (tyuatd/spectrum) 

— 19 \AAA/ 



,(3) 



X 



reduced/ 
diagrams 





o 







tree diagrams 



FIG. 10: This shows the correspondence between reduced, diagram and the tree diagram for the fifth-order spectrum as an 
example. The arrows indicate the focused vertex in each diagram. 



(CkiCksCkaCkXksK 



NL 



(P c (fc 1 )P c (fc 12 )P c (fc 4 5)P c (fc 5 ) + 59 perms.) 
-uf L (P C (^i)^c( fc i2)^c( fc 4)Pc( fc 5) + 59 perms.) 
(P c (fci)P c (fc 2 )P c (fc 3 ) ^c(^) + 4 perms.) , 



(IV.8) 



with 



NL 



l {2) 
L NL 



L {3) 
NL 



N F 



dNN a (N)Qf 3)bc (N)n b (N)n c (N)j , 
dN I N a (N)Q^ bcd (N)Q b (N)N c (N)N d (N) + 3Q a (N)Q<{ 3)bcd (N)n b (N)N c (N) 1 , 



N F 



dN^N a (N)Ql 5)bcde (N)N b (N)N c (N)N d (N)N e (N) + 6<S> a (N)Q? 3)bc (N)N b (N)N c (N) 
+AQ a (N)Q a i4)bcd (N)N b (N)N c (N)N d (N) + 12* a (N)Q a (3)bc (N)N b (N)N c (N) 

+3N a (N)Q a i3)bc (N)fl b (N)fl c (N) 1 . (IV.9) 
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Here we introduced new vectors defined by the equations 

d ~ 



^W(N) = P a b (N)n\N), (IV.10) 

^a(N) = -P b a (N)MN) - N b (N)Q b 4)cda (N)N c (N)N d (N), (IV.ll) 

_d_ 

dN 
_d_ 

dN 



^.(AO = -P b a (N)V b (N) - n b (N)Q\ 3)ca (N)N c (N), (IV.12) 
J-fl«(N) = P a b (N)fl\N) + Qf 3)bc N b (N)N c (N), (IV.13) 



(IV. 14) 

with the boundary conditions O a (A») = A ab Q b (N*), <S>a(N F ) = N[ bc N b (N F )N c (N F ), # a (iV» = and n Q (A») = 0. 

V. DISCUSSION AND CONCLUSION 

The primordial non-Gaussianity has been focused on by many authors as a new probe of the inflation dynamics. 
The deviation from the Gaussian statistics affects not only the bi-spectrum of the primordial curvature perturbations 
but also the higher order correlation functions. In general, to describe the higher order correlation functions, we need 
more parameters and more complicated calculations. Instead of resorting to the direct calculations, we developed a 
diagrammatic method, which is useful in counting the number of necessary non-linearity parameters and computing 
the higher order correlation functions for non-Gaussianity of local type. We showed that the number of parameters 
to describe the n-point correlation function is equal to the number of reduced tree diagrams with n vertices that are 
not isomorphic to each other. We also found that in the calculation of general rt-point correlation function we have 
only to solve the vector quantities which follow the same linear perturbation equation for the background field or it's 
dual pH . but with a source term and different boundary conditions. Our formalism requires the number of operations 
proportional to M even for higher order correlation functions, in contrast to the naive expectation oc J\f n , where M 
is the number of components of the inflaton field. It will be clear that our formulation is particularly powerful for 
the inflation models with many components of scalar field, including the models in which the slow-roll conditions are 
violated after the horizon crossing time. 

In this paper, we assumed that the distribution of initial perturbations of the field {£np%} is Gaussian. As a results, 
in the diagram the number of lines connected to the open circle, which corresponds to the contraction of 5<p%, is two. 
When we need to consider the effects of non-Gaussianity of Sip*, we can easily extend our formalism by adding open 
circles with appropriate numbers of the attached lines. For example, it has been well known that the leading effect of 
non-Gaussianity in 5<p® affects the three-point correlation function as [12j, [l9|, l20j 



(CkiCkaCks, 



N* a N* b N* c B abc {k u k 2 , fc 3 ) + (P c (fc 1 )P c (fc 2 ) + 2 perms.) 



5( 3 >(k 1+ k 2 +k 3 ) , (V.l) 



{StfSiplSrf) = B abc (k 1 ,k 2 ,k 3 )5^(k 1 +k 2 +\a 3 ) , (V.2) 

where denotes the non-linearity parameter given by Eq. (|IV.4j) . which has been obtained under the assumption 

that 5(p® is Gaussian. In our diagrammatic method, the first term on the right hand side of Eq. (IV. ip can be described 
by the diagram presented in Fig. QT] For this open circle o with three legs we assign the factor B abc defined in (fV2|) . 
Generalization of taking into account the higher order correlators is straightforward. Application of our formulas to 
some explicit models will be reported in the forthcoming paper. 
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FIG. 11: This diagram corresponds to the leading order correction due to the non-Gaussianity of the initial perturbation <5</?J 
to the bi-spectrum. 
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APPENDIX A: PROOF OF SUBSECTION ITTBl 

We give a proof of the statement that the functions obtained by applying the rules in Sec. Ill B | to two tree diagrams 
with n vertices that are not isomorphic to each other show different wavenumber dependence 23]. To prove it, it is 
enough to show that we can uniquely reconstruct the tree diagram with n vertices from a given function /(fei, • • • , k n ), 
which guarantees one-to-one correspondence between a diagram and a function f(ki, ■ ■ ■ ,k n ). Here each wavenumber 
ki is assigned to each vertex. By construction, the function / should be a product of (n— 1) power spectra, P, whose 
arguments arc the length of the sum of several wavenumbers taken from fcj(l < i < n). 

Let us focus on one arbitrary vertex of the would-be reconstructed diagram. We refer to this vertex as V m and the 
vector attached to this vertex as k m . We eliminate k m from the arguments of P by using the relation k\ + • ■ • + k n = 0. 
Then, the wavenumber assigned to a vertex connected to V m by a line must appear in / only once. This is because 
such a wavenumber appears only in P corresponding to the line that connects this vertex to V m . By finding all such 
wavenumbers, we recognize all the vertices that are connected to the vertex V m . By doing the same thing for each 
vertex, we completely recognize how all the vertices are mutually connected. Obviously, this fixes the shape of the 
diagram uniquely. 



APPENDIX B: EXPLICIT FORMULAS FOR DERIVATIVES OF 7V f 



As mentioned in Sec. IIII Al N F ia2 a ^ defined in Eq. can be written in terms of local quantities evaluated at 

N = Np. Here, as examples, we explicitly evaluate the coefficients N^,N[ b and N^ bc . Taking the hypersurface at 
N = Np to be a uniform Hubble one, which is equal to the uniform density slicing on super-horizon scales, we have 
the equation, 



H {y a {N F + C(N F ))) = H(<f a (N F )) 



(B.l) 



The Hubble parameter H is given by Eq. pi.5[) . In our previous paper [ll[, solving Eq. (|B.1[) with respect to ((Np) 
up to the second order, we have obtained 



N 



N F h = - 

1 * ab 



H b (<p)F»(tp) 



H c {y)F c {v) 



(0) 

v= tp(N F ) 



<°>, 

<p= V (JV» 



where 



H c P c d F d ) « , 



U ab = H ab + 2 (H c P c a + F c H ca ) N£ + (F c H cd F 
with H a = dH/dip a , FL a b = d 2 H/d(p a dip b . Solving Eq. (|B.1[) up to the third order, we also obtain 



N: 



W abc (ip) 



be 



H d (<p)F*(<p) 



(0) 

v =ip(N F ) 



(B.2) 
(B.3) 

(B.4) 

(B.5) 
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where 



W a 



be 



H a bc 

+3 



H d (Q% )ef F e + P d e P e f ) F f + H def F d F e F f + 3F d H de P e f F f 



2F d H de P e a + (H ade F d + H ad P d e ) F e + H d [Qf 3)ea F e + P d e P e a 

-3 (2H ad P d b + F d H dab + H d Q d 3)ab ) N? 

-3 (F d H de F e + H d P d e F e ) N[N b F c + 3 [F d H da + H d P d a ) N b F c 



N^N b F N[ 



N b F Nf, 



(B, 



with H abc = d 3 H / d(p a dip b d(p c . Note that U a b and W a bc are symmetric with respect to the indices. As we mentioned 
earlier, here we define the phase space variables as ip{ = ft and ft 2 — ft ■ 
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